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Abstract 

Background: Recovering the network topology and associated kinetic parameter values from time-series data are 
central topics in systems biology. Nevertheless, methods that simultaneously do both are few and lack generality. 

Results: Here, we present a rigorous approach for simultaneously estimating the parameters and regulatory 
topology of biochemical networks from time-series data. The parameter estimation task is formulated as a mixed- 
integer dynamic optimization problem with: (i) binary variables, used to model the existence of regulatory interac- 
tions and kinetic effects of metabolites in the network processes; and (ii) continuous variables, denoting metabolites 
concentrations and kinetic parameters values. The approach simultaneously optimizes the Akaike criterion, which 
captures the trade-off between complexity (measured by the number of parameters), and accuracy of the fitting. 
This simultaneous optimization mitigates a possible overfitting that could result from addition of spurious regulatory 
interactions. 

Conclusion: The capabilities of our approach were tested in one benchmark problem. Our algorithm is able to 
identify a set of plausible network topologies with their associated parameters. 
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Biochemical networks 



Background 

Mathematical models of biochemical systems are be- 
coming essential in systems biology to complement and 
extract information from time series. This information 
can be of two types. On the one hand, if the structure of 
the molecular circuit that executes the process of inter- 
est is known, models can be used to infer the numerical 
parameters that govern the dynamics of the system 
[1-4]. On the other, models can be used to infer the 
structure of the system from time series data (see for 
example [5-7]). 

In either case, to obtain a useful model, we face different 
challenges: (i) defining the system's mass flow structure 
{stoichiometry), (ii), deciding the appropriate mathemat- 
ical representation (kinetics), (iii) estimating the parame- 
ters that make the model response consistent with 
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experimental data (parameter estimation), and (iv) infer- 
ring the system's regulatory structure. In addition, once 
the model is well defined, it should be able to predict 
systemic responses under yet untested experimental 
conditions (model validation). 

The four challenges described in the previous para- 
graph are often addressed in independent steps. Current 
solutions to the first challenge are generally based on 
compiling information about the system and using that 
information to create the stoichiometric matrix for the 
system one wants to analyze (see for instance [8]). To 
solve the second challenge we need to define kinetic 
functions that describe the dynamic behavior of the 
dependent variables of the system. If the kinetic func- 
tions are unknown, approximate formalisms that have a 
solid theoretical support can be used to describe the dy- 
namic behavior of the system within a given accuracy 
[9,10]. The third challenge is typically formulated as an 
optimization problem that minimizes the sum of 
squared residuals between the measured and simulated 
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data (see a review of methods in [1]). The type of 
optimization problem being faced and the technical 
challenges to be solved depends upon the biological 
model of choice, upon the experimental data available, 
and upon the specific mathematical formalism used 
[11,12]. In many practical applications, the target bio- 
logical system is described through nonlinear ordinary 
differential equations (ODEs). Hence, the parameter esti- 
mation task gives rise to dynamic optimization problems 
that are hard to solve. The fourth challenge could in 
principle be addressed in the same way as the first. How- 
ever, despite the enormous amount of biological infor- 
mation available in public databases, regulatory signals 
are, in general, poorly understood and hardly ever prop- 
erly characterized in vivo. Regulatory signals appear in a 
model as parameters accounting for the influence that 
metabolites others than the substrates of a reaction have 
on its velocity. Hence, parameter fitting can also be used 
to address the fourth challenge. However, the over- 
whelming majority of parameter estimation methods as- 
sumes a given structure and considers a fix regulatory 
scheme (see a review in [1]). This simplification is moti- 
vated by the difficulty in identifying regulatory effects, a 
task for which a myriad of alternative kinetic models 
must be explored [7,13-15]. 

Traditional methods for the selection of biological sys- 
tems have mostly applied regression or chi-squared- 
based criteria (rather than information-theoretic fit 
criteria) [16]. However, information-theoretic criteria 
such as the Akaike's Information Criterion (AIC) [17] or 
the Bayesian Information Criterion (BIC) [18], are now 
perceived as important measures to assess quality of fit. 
AIC is often preferred over BIC becaue it has a more 
immediate connection to the theory of information [19]. 
AIC captures the trade-off between the complexity 
(measured by the number of parameters), and accuracy 
of the fitting. Smaller AIC values imply a better approxi- 
mation to the model sought. 

In this work we propose a strategy to simultaneously 
address the four challenges described above that relies 
on the use of mixed-integer dynamic optimization 
(MIDO) methods. Our approach adopts a structured 
mathematical framework to represent the kinetics of the 
processes that is flexible enough to reproduce a set of 
plausible network topologies (by implementing slight 
modifications on a basic model formulation). The power- 
law [20] and the saturable and cooperative formalisms are 
examples of such general kinetic representations [9]. 
Based on this type of general kinetic modeling framework, 
we develop our systematic parameter estimation method 
that provides as output a set of potential reaction and 
regulatory topologies for the target network along with 
the associated model parameters. We illustrate the 
capabilities of our approach using the GMA kinetic 



representation, a canonical model structure that uses the 
power-law kinetic formalism [21,22]. 

Results and discussion 

As a proof-of-concept, we have tested the capabilities of 
our approach through its application to a case study 
taken from Voit and Almeida [23]. The system consid- 
ered is a four-constituent pathway branched with six vel- 
ocities and two regulatory signals. Xj is generated from 
X 0 , and its production is inhibited by X 3 which is pro- 
duced from Xj via intermediate X 2 . Xj yields also X 4 , 
which promotes the degradation of X 3 (see Figure 1). 

Parameter estimation when the regulatory structure 
is known 

We shall first show that the proposed method is capable 
of appropriately identifying the model parameters using 
dynamic data when the regulatory structure is known. 
This is the classical parameter estimation problem that 
is solved in many applications. To this end, we first pro- 
duce dynamic data without error from the reference sys- 
tem using a specific set of parameter values. Then, this 
in silico data is labeled as experimental and we use the 
proposed method to estimate the model parameters. We 
define a dynamic optimization model that contains a set 
of dynamic differential equations describing the system's 
kinetics. This dynamic model is reformulated into a non- 
linear program (NLP) using orthogonal collocation on 
finite elements. This NLP does not contain binary vari- 
ables because we assume that the regulatory signals are 
known. The aforementioned NLP was implemented in 
GAMS 23.7.3 and calculated with CONOPT 3.15A on a 
PC/AMD Athlon at 2.99 Ghz using a single core. The 
NLP features 302 variables and 285 constraints, and was 
solved in 2.3 CPU seconds. As expected, we obtain 
estimated parameters values that are very close to the 
original ones (see Table 1), and a least square error of 
1.45 x 10' 6 . 

Non-linear kinetic models, like the GMA representa- 
tion, have a certain degree of plasticity that allows differ- 
ent parameter sets to fit the same data. Clear parameter 




Figure 1 Reference system taken from Voit and Almeida [23] 

(default parameters are shown in Table 1 ). 
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Table 1 Original and predicted parameters values 



Parameter 


Original parameters 


Proposed algorithm 




-0.8 


-0.7999 


h\ 


0.5 


0.4996 


hi 


0.75 


0.7494 


Ui 


0.5 


0.5006 


fs3 


0.5 


0.4996 


Ua 
'M 


0.2 


0.1996 




0.8 


0.8010 


Yi 


12 


12.000 


Y2 


8 


8.0031 


Y3 


3 


3.0034 


T4 


2 


1 .9965 


T5 


5 


5.0014 


T6 


6 


5.9967 



Data is error free (one experiment with only one observation by time point). 



trends are obtained by fixing a given parameter and fit- 
ting the remaining ones. As an example, Figure 2 shows 
the results of fixing f 32 at different values and fitting the 
other parameters. All the points in the figure lead to re- 
siduals below 5.88 x 10~ 4 , indicating that it is possible to 
obtain good fits with different parameter sets. Similar 
patterns are obtained if we choose to fix any other par- 
ameter of the set. 

As observed, the model is rather flexible, as there are 
many combinations of parameters values leading to very 
low residuals and essentially the same fit to the data. In 
practical terms, this means that given an experiment and 
an estimation procedure, we could obtain different par- 
ameter sets that closely reproduce the experimental 
measurements, but that differ from the actual values 
with which the dynamic profile was generated in silico. 
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Figure 2 Values of the fitted parameters for different values of 
f32. Each point was generated by fixing f32 and solving the NLP 
free of error. 



Thus, estimated parameter values don't help comparing 
the obtained fit with the reference model. In practice, 
the residual error and the resulting time profiles should 
be used to assess the fit. 

We will now consider the effect of noisy data on fitting 
the model, as such noise plays a key role in evaluating 
any proposed method for identifying the regulatory 
structure of a network. To explore the influence of ran- 
dom experimental uncertainty, we generated 100 dy- 
namic profiles from the reference model by introducing 
statistical noise. For this, we applied Monte Carlo sam- 
pling assuming that every data point follows a normal 
distribution with standard deviation values of 0.5, 1, 5 
and 10% of the actual nominal value. For comparison 
purposes, we use the same perturbation experiment as 
in the previous example. Table 2 shows the parameter 
values and the associated residuals obtained for four of 
the samples generated, while Figure 3 depicts the pro- 
files associated with a standard deviation of 10%. We 
can appreciate that despite the different parameter 
values, the various fitted models lead to similar residuals. 
Note that although the regulatory structure is fixed, we 
obtain parameter values representing either positive or 
negative regulatory effects (f S4 ) of X 4 on v s . This is a 
consequence of the "experimental error" introduced in 
the noisy data. That error may force the estimation pro- 
cedure to an optimum involving a set of parameter 
values that may be different from the set that generates 
the noiseless data. In addition, as seen above, different 
parameters sets can be used to produce similar time 
courses. This means that there are coupled parameters 
in the system, which may also contribute for the estima- 
tion of regulatory interactions with reversed signals. 

In general, even in simple cases as the one considered 
here, it will be difficult to obtain a consistent estimation 
from a single time-series. Identifying the parameter set 
that is more likely to be the correct one requires simul- 
taneous fitting to additional time-series, resulting from 
more than one set of experiments. By doing so, we will 
constraint further the admissible parameter sets (see 
[24]). In Table 3, we show the results of fitting three dif- 
ferent experiments with experimental error. Each experi- 
ment corresponds to an alternative perturbation on the 
initial concentration of metabolite X 3 (0.2, 1.2, and 2.2). 
These perturbations force the system to move across dif- 
ferent dynamic regimes, producing additional informa- 
tion that helps in the identification of appropriate 
parameter values. As observed, the estimated parameters 
are more consistent over the various experiments. They 
are also closer to the actual parameter set selected for 
generating the data. Note, however, that it is still pos- 
sible to find solutions involving alternative regulatory 
topologies with good fit to data (f 54 acting as an inhibitor 
in Profile 2). 
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Table 2 Parameters values with noisy data (one experiment) 



10% 





Profile 1 


Profile 2 


Profile 3 


Profile 4 


f» 


-0.14 


-0.27 


-0.84 


-0.79 


fx 


0.26 


0.47 


0.4 


0.29 


hi 


0.44 


1 


0.64 


0.41 




0.04 


0 


0.9 


1 


fa 


0 


0.26 


0.42 


0.12 




-0.06 


0.04 


0.1 


-0.12 


'64 


0.13 


0.07 


1 


1 


Residua 


1.88 


1.67 


1.68 


2.29 


5% 




Profile 5 


Profile 6 


Profile 7 


Profile 8 


f"l3 


-0.282 


-0.532 


-0.631 


-0.893 




0.56 


0.618 


0.306 


0.6 


hi 


1 


1 


0.436 


1 




0 


0.092 


0.761 


0.742 


fa 


0.368 


0.639 


0.273 


0.298 




0.127 


0.244 


0.021 


0.279 


f(A 


0.064 


0.158 


I 


1 


Residua 


0.41 28 


0.4203 


0.5706 


0.4482 


1% 




Profile 9 


Profile 10 


Profile 1 1 


Profile 12 




-0.881 


-0.427 


-0.859 


-0.71 


hi 


0.571 


0.523 


0.5 


0.414 


hi 


0.885 


0.809 


0.758 


0.608 


U, 


0.587 


0.078 


0.661 


0.656 


fa 


0.479 


0.467 


0.507 


0.402 


fs4 


0.2 


0.176 


0.197 


0.136 


f(A 


1 


0.162 


I 


1 


Residua 


0.0207 


0.0163 


0.0167 


0.0227 


0.5% 




Profile 13 


Profile 14 


Profile 15 


Profile 16 


fa 


-0.845 


-0.744 


-0.843 


-0.765 


fx 


0.535 


0.472 


0.496 


0.453 


hi 


0.816 


0.714 


0.749 


0.673 




0.556 


0.492 


0.647 


0.643 


fa 


0.492 


0.439 


0.497 


0.443 


fa* 


0.201 


0.167 


0.196 


0.164 


fa 


0.916 


0.816 


1 


1 


Residua 


0.0052 


0.0041 


0.0042 


0.0057 



We solved a total of 100 problems, each corresponding to a different 
replication, generated randomly see Additional file 1: Table SI). The table 
shows the 16 cases for which the residual error is low. 



Identifying the regulatory structure 
Performance using error free data 

After testing the capabilities of the method when the 
structure is known, we studied its ability to identify the 




Time 



Figure 3 Adjusted profiles for four different noisy data sets (i.e. 
one experimental condition and four replications) with a standard 
deviation of 10%. 

V / 

regulatory topology of the model. To this end, we ex- 
plore the performance of the method using one experi- 
ment with low experimental error (i.e., assuming that 
the data follow normal distributions with a standard de- 
viation of 0.5%). Larger errors result in a wider set of 
alternative structures and for simplicity's sake we shall 
not discuss them here. 

In order to simplify the search, we fix a maximum of 
two metabolites (the substrate of the reaction, which is 
given by the stoichiometric information, and one pos- 
sible additional modifier, which is not a priori character- 
ized) as potential variables affecting each velocity. 

We note that it is typical to have some a priori know- 
ledge about the biological system one is interested in. 
The complexity of the regulatory interactions in the 
identification problem is reduced if such knowledge can 
be used to constrain further both, the number of poten- 
tial regulatory signals in the model and their signs (posi- 
tive, negative). In such cases, we can introduce specific 



Table 3 Parameter values obtained from simulated noisy 
data (with noisy data (three experiments)) 





Profile 1 


Profile 2 


Profile 3 


Profile 4 


fa 


-0.67 


-0.64 


-0.62 


-0.92 


fa 


0.33 


0.9 


0.49 


0.69 


hi 


0.42 


I 


0.73 


1 


fa 


0.64 


0 


0.38 


0.26 


fa 


0.49 


0.66 


0.3 


0.4 


fa 


0.05 


-0.95 


0.22 


0.34 


fa 


I 


1 


0.53 


0.58 


Residual 


6.96 


7.10 


5.39 


4.89 



We solved a total of 100 problems, generated randomly. See Additional file 1: 
Table S2. 
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constraints for the relevant parameters to be fitted. For 
example, in our case kinetic-order corresponding to the 
substrates of a reaction must be positive. 

The MINLP model that simultaneously fits the param- 
eters and infers probable regulatory interactions was im- 
plemented in GAMS 23.7.3 and solved with the solver 
SBB in the same computer as before. The model has 72 
binary variables, 391 continuous variables and 414 equa- 
tions. The solution time was in the order of few minutes 
for each simulation. 

Our algorithm identifies a set of compatible systems, 
since the model has enough flexibility to play with the 
regulatory structure as well as the kinetic parameters 
when minimizing the residuals. The method identifies 



Residual = 0.00223 
! C2)sr Xo Xo ^"^3 

Residual = 0.00283 



© 
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Residual = 0.00316 



©. 
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Figure 4 The proposed method identifies different regulatory 
topologies that essentially produce the same output. We show 
here the associated profiles corresponding to three regulatory 
structures with lowest residual values obtained by analyzing data 
from a single experiment with one replicate (see parameters values 
and residuals in Additional file 1: Table S3. 



topologies that are quite close and that show very small 
residuals, but it is unable to uniquely identify the ori- 
ginal topology (see Additional file 1: Table S3 for a list of 
topologies generated and their associated kinetic param- 
eters and residuals). As an example, in Figure 4, we 
compare three completely different regulatory structures 
that produce almost indistinguishable results and similar 
fitting to the actual dynamics, leading to residual values 
of 0.00223, 0.00283 and 0.00316 (Figure 5). 

As before, one strategy for increasing the possibility of 
correctly identifying the "true" regulatory structure is to 
use additional time-series data of the same system under 
different sets of initial conditions. To this end, we chan- 
ged the initial concentrations of X 3 (0.2, 1.2, and 2.2). 
The MINLP model was again implemented in GAMS 
and solved with SBB in the same computer. In this case, 
the MINLP features 72 binary variables, 967 continuous 
variables and 980 equations. The solution time was in 
the order of few minutes for each simulation. 

In Figure 6 we show the dynamic profiles associated 
with three different topologies identified by the MINLP. A 
complete list of network topologies and associated kinetic 
parameters and residuals is provided as (Additional file 1: 
Table S4). With three time series, the method identifies 
not only the actual topology, but also several structures 
that contain the original one (i.e., topologies that account 
for all the actual regulatory effects plus other signals that 
were not present originally). Again, we obtained slightly 
different parameter sets in each case, since the model 
flexibility is rather large. 

Additional remarks 

The use of MIDO techniques combined with orthogonal 
collocation allows posing the parameter estimation task 
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Time 

Figure 5 Dynamic responses corresponding to the three 
different topologies of Figure 4. Parameter values are indicated 
on Additional file 1: Table S3. 
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Time 



Figure 6 The Profiles generated from three different topologies 
and three experiments with one replication each. The experiments 
are generated from the base case by applying different perturbations 
in the initial concentration of X 3 . Details on the topology and 
associated parameters are provided on Additional file 1: Table S4. 



as an algebraic optimization problem that can be efficiently 
solved using standard MINLP algorithms. Orthogonal col- 
location shows some appealing properties (see [25]), but 
has the drawback of increasing the model size because it 
adds auxiliary variables and equations that increase the 
problem complexity. Our MIDO approach, however, can 
be solved by any MIDO algorithm, and it is not re- 
stricted to the use of orthogonal collocation and MINLP 
reformulations. 

A key point in our method is the selection of an ap- 
propriate starting point to initialize the MINLP algo- 
rithm. Standard MINLP algorithms typically solve an 
initial NLP where the binary variables are relaxed. If this 
NLP does not converge, the entire algorithm might fail. 
An initialization strategy that works well in practice is to 
integrate first the original kinetic model for some par- 
ameter values, and then use the dynamic profiles gener- 
ated in silico to provide a starting point for the NLP 
solver. Another method consists of solving an auxiliary 
model where we relax some constraints through the 
addition of slack variables, and then minimize the sum- 
mation of the slacks in order to obtain an initial feasible 
point. With this relaxed model, we can identify a feasible 
(but not necessarily optimal) solution for the initial NLP. 

Even if the MINLP model converges, there is still the 
issue of getting trapped in local optima during the 
search. To avoid this, we can run the optimization algo- 
rithm from different starting points generated randomly. 
This strategy does not guarantee convergence to the glo- 
bal optimum, but tends to produce high quality solu- 
tions in short CPU times. In contrast, deterministic 
global optimization methods provide a rigorous interval 



within which the optimum should fall, but tend to lead 
to large CPU times (see [26,27]). 

In our case, we initialize the NLPs by solving a set of 
relaxed problems from different starting points and then 
pass these results to the first NLP solver. This approach 
provides feasible points from which the model converges 
to solutions with low residuals. 

In general, due to the nonconvex nature of the refor- 
mulated MINLP, the nonlinear branch and bound imple- 
mented in SBB outperforms the outer-approximation 
used by DICOPT. This is because the supporting hyper- 
planes defined in the master MILP solved by DICOPT 
may chop-off feasible solutions due to the noconvex 
nature of some nonlinear inequalities. 

We note that nonlinear models are hard to handle, 
and even more so when they contain binary variables. 
Standard NLP solvers can solve problems containing up 
to hundreds of thousands of variables and constraints. 
On the other hand, the computational burden of MIDO 
(and MINLP) models is rather sensitive to the number 
of binary variables. For the type of problems we are deal- 
ing with, it is difficult to provide a bound on the number 
of binaries above which the algorithm might fail. In prac- 
tice, however, we found that this approach efficiently for 
less than one hundred binaries (around 30 parameters). 

From a practical viewpoint, we face the challenging 
problem of discriminating between compatible regula- 
tory structures for a given data set. On a worst case sce- 
nario, our method provides a ranked set of alternative 
regulatory topologies that can be tested and validated 
experimentally. If appropriate additional time-series data 
are available, the set of admissible solutions for testing 
can be further constrained and reduced. Our method 
finds a set of alternatives that are consistent with the dy- 
namic data available and that can be further refined 
using additional information and expert knowledge on 
the system, (i.e., complementary biological information). 
For instance, kinetic-orders that correspond to substrates 
of a reaction may be safely restricted to be positive. Simi- 
larly, if we are fairly sure that a given metabolite does not 
participate in a reaction, its kinetic-order should be fixed 
to zero. 

Our method can also be used to explore hypotheses 
about the regulatory structure of a system. For instance, 
we can force some parameters to take negative values, 
thereby representing inhibition effects, and then perform 
the optimization so as to determine if the fitting is good 
enough. Furthermore, we can follow the same procedure 
in order to identify regulatory effects that are consistent 
with this hypothesis. 

In addition, we note that our approach can be easily 
adapted in order to work with other model selection cri- 
teria besides AIC. We remark, however, that the assess- 
ment of different selection criteria would deserve a 
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comprehensive study that is beyond the scope of this 
work. 

The simple examples presented in this paper show that 
estimating parameters in dynamic kinetic models is far 
from being an easy job and that models based on the 
power-law formalism facilitate the estimation task. Al- 
though this formalism is suitable for a wide variety of 
problems, one may argue that it may present some limi- 
tations. As an alternative, we can use extensions of this 
framework such as the Saturable and Cooperative for- 
malism [9], which takes into account saturation effects. 
In both cases, a key point is the possibility of using a ca- 
nonical mathematical formalism that facilitates the auto- 
matic search of alternative regulatory patterns. The 
method described here would be applicable to such 
models via recasting of the Saturating and Cooperative 
formalism into a power law [28]. 

Conclusions 

In this work we have proposed a rigorous approach 
based on mathematical programming for the simultan- 
eous identification of the regulatory signals and estima- 
tion of the kinetic parameters of models of biochemical 
networks. Our approach is based on the use of mixed- 
integer dynamic optimization (MIDO) models that 
minimize the Akaike criterion, and that can be solved by 
standard optimization algorithms. Particularly, we solve 
this MIDO by reformulating it as a mixed-integer non- 
linear program (MINLP) using orthogonal collocation 
on finite elements, which makes it possible to apply 
standard MINLP solution algorithms in an iterative fash- 
ion in order to identify a set of plausible network topolo- 
gies and associated kinetic parameters. 

It is noteworthy that the difficult task of parameter es- 
timation in nonlinear models becomes really compli- 
cated as the size of the models increases. Therefore, 
such estimation typically requires customized solution 
procedures. One key point is to use the appropriate ini- 
tial conditions to ensure convergence of the calculations. 

The proposed method can contribute to fill the lack of 
information on the regulatory signals that are in play in 
a given metabolic scenario. Although we cannot deal 
with genome-wide models, we have shown that dynamic 
profiles can be processed to provide clear hypothesis on 
the underlying regulatory structure. This is an important 
step towards completing essential information on differ- 
ent metabolic processes that are poorly understood. 

Methods 

The problem we address here is to infer the regulatory 
structure of a metabolic system, given a known structure 
for the reaction network (stoichiometry) and experimen- 
tal time series for the dynamic behavior of that system. 
To address this question, and to explore the practical 



problems associated, we consider the following general 
representation of a biochemical network: 

p 

*i = ^2fi,r V r i= l,-,n (1) 
r=l 

where X t denotes the concentration of metabolite i, fi iir 
is the stoichiometric coefficient of metabolite i in 
process r, which indicates the number of molecules of type 
i produced or destroyed by process r, and v r is the rate 
function of this process. In general, v r is represented as: 

v r (Xi, X n+m . 6) (2) 

There are two critical issues in defining this model. 
One is the selection of an appropriate mathematical rep- 
resentation for v„ which may be a function of an arbi- 
trary number of variables (substrates, products, and 
modifiers). In most cases the mechanism for each 
process are unknown and choosing a specific mechanis- 
tic rate law, such as a Michaelis-Menten rate law, be- 
comes an act of faith. The other issue is the problem of 
identifying the regulatory structure of the system. 

The most straightforward and theoretically well sup- 
ported solution to both issues is the use of an approxi- 
mate formalism based on a standard mathematical 
representation [10]. By adopting such a kinetic represen- 
tation, identifying the regulatory structure of the system 
becomes synonymous to determining the set of values 6 
for the model parameters that better fit the available 
data. Hence, without losing generality, and as a first step 
towards a more complex framework, we will consider 
the case where the rates are modeled using a power-law 
formalism. Note, however, that our approach could be 
easily extended in order to accommodate any other 
structured kinetic formalism. 

Power-law models 

Using the power-law representation, the rate v r is ex- 
pressed as follows: 

n+m 

Vr= Y r J{Xj r = 1,.., p (3) 

//I 

where y r is an apparent rate constant for reaction r, and 
f n j is the kinetic order of metabolite / in that process. 
Note that this equation accounts for the effect of n + m 
metabolites (n dependent and m independent) on each 
reaction. 

The advantage of this representation is that the same 
functional form represents all the rates. The reaction 
structure of the system will constrain the range of ad- 
missible values for some of the parameters. For example, 
all y and /parameters for the substrates and catalysts of 
the reactions are by definition larger than zero. In 
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addition, the values of the / parameters for all metabo- 
lites that are not directly involved in a given process are 
zero in the rate that describes the process. 

By adopting such a kinetic representation, we can pose 
the problem of identifying the regulatory signals in a 
very compact mathematical form. If Xj is a modifier of 
v n then the corresponding kinetic order f r j will be differ- 
ent from zero (positive if it is an activator, and negative 
if it is an inhibitor). By substituting (3) into equation (1), 
we get what is known as a Generalized Mass-Action 
(GMA) model. 



Srj 



(4) 



Note that the power-law formalism accounts for both 
the stoichiometry of the system (the network structure), 
and the reaction and regulatory structures (kinetic or- 
ders) using a single systematic nonlinear representation. 
This property is very important for defining a systematic 
way of exploring alternative regulatory signals. We will 
make use of this general and compact formalism in the 
derivation of the equations for the parameter estimation 
model. 

Parameter estimation in a GMA model 

Given a set of experimental observations (i.e., time 
courses for the metabolites), our goal is to find the 
values of the apparent constants and kinetic orders that 
minimize the sum of least squared errors between the ex- 
perimental data and the predicted dynamic profiles. This 
problem can be expressed in compact form as follows: 



min 



s.t. 



Ii=l 1=1 

p 

x i = ' = 1 



r-1 
n+m 



Vr = JvIPv r = 1 '->P 



i-i 



Xi(t 0 )=X oi i=l,...,n;t e[f 0 ,f/] 
Xff = X t (t u ) i = !,...,«;« =1,...,* 



(5) 



where X t represents the state variables (i.e., metabolite 
concentrations), X oi their initial conditions, X^f denotes 
the experimental observations, and X%% are the values 
calculated by the dynamic model (i.e., model predic- 
tions), i is the index for the set of state variables whose 
derivatives explicitly appear in the model, y r and f r , are 
the parameters to be estimated, and t u , is the time asso- 
ciated with experimental point u belonging to the set U 
of observations, k is the total number of experimental 



data points and n is the number of time dependent 
variables. 

Conventional parameter estimation approaches seek 
parameter values that minimize the approximation error 
assuming a given regulatory scheme (i.e., fixing some f n j 
to zero beforehand according to the aprioristic biochem- 
ical knowledge of the system). While this assumption 
simplifies the calculations, it can lead to poor approxi- 
mations and hamper at the same time the discovery of 
new regulatory loops. In this work we introduce a rigor- 
ous and systematic parameter estimation and network 
identification method that makes no assumption regard- 
ing the regulatory network topology. 

To model the existence of a regulatory interaction, we 
make use of the following disjunction: 



V 



V 



Y+. 

rj 



j= 1, ...,«, 

r = l,...,p 
(6) 



In which Y r J,Y r ,i and Y r f are Boolean variables that 
are true if parameter f r i is negative, zero or positive, re- 
spectively, and false otherwise, e is a very small param- 
eter. Note that only one term of the disjunction can be 
active (i.e., exclusive disjunction), while the others must 
be false. Hence, if Y r i is true, metabolite i takes no part 
in velocity r. Conversely, if this metabolite has an influ- 
ence on r, then Y r j is false and either Y r J or Y r t will be 
active. This disjunction can be translated into standard 
algebraic equations using either the big-M or convex- 
hull reformulations [29]. By applying the former, we get: 

frj * - E + M { l -y~r,i ) /=!,•», n,r=l,...,p 

-s-M{l-y rj j <f rj ze + M[l-y rj j j=l n,r=l,...,p 

f rj < e + M(l-y+j ) ; = 1,..., n,r=l, ...,p 

yrj +y r ,j +ytj =1 ; = 1 >-> n,r = i,...,p 
y;,i +y r ,j +fi e {0,1} 

(7) 

where Boolean variables Y have been replaced by auxil- 
iary binary variables y. In these equations, M is a suffi- 
ciently large parameter whose value must be carefully 
set according to the bounds defined for the kinetic 
parameters. 

A key issue in our approach is how to avoid overfit- 
ting. To this end, we make use of the Akaike criterion, 
which captures the trade-off between the number of kin- 
etic parameters contained in the model and its ability to 
accurately reproduce the experimental data. If we as- 
sume that the error of the observations follows a normal 
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distribution, the Akaike criterion takes the following 
mathematical form [17]: 



/ k " f s2\ 
I \ " \ " / v exp v mod \ 



AIC= /clog 



V 



2 EE(^,/ + C 

;=1 r=l 



(8) 



Where AIC denotes the value of the Akaike criterion 
and C is a constant value that does not affect the 
optimization. The parameter estimation problem can be 
finally posed in mathematical terms using the following 
MIDO (mixed-integer dynamic optimization) formulation: 



(M) min k log 



\ 



J 



S.t. 



• 2 EEU ••'.) 

p /=! r=l 

X) = j2Pi,r V r » =!,•»,» 
r-1 

vr = r r il x j r=1 >->P 

Xi(t 0 )=X 0i i=l,...,n;t e[t Q ,t f ] 
Xi M = Xi(t u ) i=l,...,n;u =l,...,/c 
frj s - £ + M { l -y~rj ) /=1,-, n,r=l,...,p 

-e-M(l-y r ^<f rJ <s + M(l-y r ^ ;=1,..., n,r = X,..,,p 
f rJ < e + M\l-y+jj 7 = 1,..., «,r=l,...,p 
J^V + =1 / = 1j-) n,r=\,...,p 

y~ r ,j +y rJ +ytj e i 0 ' 1 ) 

(9) 

There are different solution methods to solve this 
MIDO (see [25]). Without loss of generality, we propose 
here to reformulate this problem into an equivalent alge- 
braic MINLP (mixed-integer nonlinear program) using 
orthogonal collocation on finite elements. This allows 
exploiting the rich optimization theory and software ap- 
plications available for MINLP in the solution of the 
MIDO. Note that the reformulated MINLP might be 
nonconvex. This will give rise to multimodality (i.e., ex- 
istence of multiple local optima), preventing standard 



gradient-based solvers from identifying the global 
optimum. Deterministic global optimization methods 
could be applied to solve the MINLP, but they might 
lead to large CPU times given the size and complexity of 
a standard dynamic problem of this type. Details on the 
application of deterministic global optimization methods 
to parameter estimation problems of small/ medium size 
can be found elsewhere [30,31]. For the reasons given 
above, in this work we will solve the reformulated 
MINLP using local optimizers. 

One important feature of our approach is that rather 
than calculating a single optimal solution, it identifies a 
set of plausible regulatory topologies by solving the 
model iteratively. That is, the model is first solved to 
identify a potential regulatory configuration represented 
by a binary solution (i.e., set of values of the binary vari- 
ables). The model is then calculated again but this time 
adding the following integer cut, which excludes solu- 
tions identified so far in previous iterations from the 
search space: 



E */+ E yrf+ E at 

(rj)eONE^ (r,/)eOJV£ it (rj)eONE+ 

- E ^ lt - E y r /- E *tS 

(r,j)^ONE u (rj)eONEi t (rj)^ONE+ 

< | ONE u + ONE it + ON |- 1 

ONE u = {(rj^yrt' 1 = 1 in the solution obtained 
in the iteration it } 

ONE it = {{r,j)\y r j lt = 1 in the solution obtained 
in the iteration it } 

ONE% = {{r,j)\y+f = 1 in the solution obtained 
in the iteration it } 

ZERO it = {(r,f) \y~^ = 0 in the solution obtained 

in the iteration it } 

ZEROu = {( r J)\y r j' t = 0 in the solution obtained 

in the iteration it } 

ZEROl = {(r,f)\yf, u = 0 in the solution obtained 

in the iteration it } 



(10) 



Where ONE it and ZEROu represent the sets of binary 
variables that take a value of one and zero, respectively, 
in iteration it of the algorithm. After adding the integer 
cut, the model is solved again to produce a new regula- 
tory topology, and this procedure is repeated iteratively 
until a desired number of configurations is generated. 
Hence, the algorithm produces as output a set of poten- 
tial network configurations (encoded in the values of the 
binary solutions) rather than a single topology. Note that 
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these regulatory topologies show a descendant value of 
the Akaike performance criterion. 
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